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ABSTRACT 

Estimates of the mass distribution and dark-matter (DM) content of dwarf spheroidal 
galaxies (dSphs) are usually derived under the assumption that the effect of the tidal 
held of the host galaxy is negligible over the radial extent probed by kinematic data¬ 
sets. We assess the implications of this assumption in the specihc case of the Fornax 
dSph by means of Wbody simulations of a satellite orbiting around the Milky Way. We 
consider observationally-motivated orbits and we tailor the initial distributions of the 
satellite’s stars and DM to match, at the end of the simulations, the observed structure 
and kinematics of Fornax. In all our simulations the present-day observable properties 
of Fornax are not signihcantly inhuenced by tidal effects. The DM component is altered 
by the interaction with the Galactic held (up to 20% of the DM mass within 1.6 kpc is 
lost), but the structure and kinematics of the stellar component are only mildly affected 
even in the more eccentric orbit (more than 99% of the stellar particles remain bound 
to the dwarf). In the simulations that successfully reproduce Fornax’s observables, the 
dark-to-luminous mass ratio within 1.6 kpc is in the range 5 — 6, and up to 16 — 18 if 
measured within 3 kpc. 

Key words: dark matter - galaxies: dwarf - galaxies: individual: Fornax - galaxies: 
kinematics and dynamics - galaxies: structure 


1 INTRODUCTION 

Over the years, the manyfold growth in the size of spectro¬ 
scopic samples of individual stars has allowed the determina¬ 
tion of accurate and spatially extended line-of-sight (LOS) 
velocity-dispersion proHles for each of the Milky Way (MW) 
“classical” dwarf spheroidal galaxies (dSphs^; e.g. Kleyna et 
al. 2001; Koch et al. 2007; Walker et al. 2007; Battaglia et 
al. 2008; Munoz et al. 2006; Battaglia et al. 2011; Walker 
2013) as well as of higher moments of the LOS velocity dis¬ 
tribution (e.g. Lokas 2009; Breddels et al. 2013; Breddels & 
Helmi 2013, hereafter BH13; Richardson & Fairbairn 2014). 
Equilibrium dynamical modeling of these data-sets leads to 
the finding that MW dSphs exhibit some of the most ex¬ 
treme dynamical mass-to-light (M/L) ratios known to-date, 
up to several 100s (M/L)©. This has triggered a large body 
of works attempting to assess the dark matter (DM) con¬ 
tent and distribution of these galaxies (for recent reviews 
see Battaglia, Helmi & Breddels 2013; Strigari 2013; Walker 
2013), which is an important piece of information to un- 
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^ We refer to the dSphs known prior to the SDSS as “classical” 
dSphs, while those discovered afterwards as “ultra-faint dwarfs”. 


derstand the nature of DM particles (e.g. Ackermann et al. 
2014) and to test the predictions of cosmological models of 
galaxy formation, for instance regarding the mass function 
of luminous satellites around MW-like systems (e.g. Klypin 
et al. 1999; Moore et al. 1999) and the effects of baryons 
on the DM distribution (e.g. Navarro et al. 1996; Nipoti & 
Binney 2015, Zolotov et al. 2012). 

However, it has always been a matter of debate whether 
the tidal interaction experienced by these small galaxies dur¬ 
ing their orbital evolution around the MW may have im¬ 
pacted their dynamical status to the extent of significantly 
altering estimates of their current DM content and distribu¬ 
tion. Observational evidence of on-going tidal disturbance 
in some M31 satellite galaxies is provided by the clear pres¬ 
ence of tidal tails or strong structural distortions (e.g. Choi, 
Guhathakurta, & Johnston 2002; McConnachie et al. 2010; 
Crnojevic et al. 2014). Around the MW, beside the spectac¬ 
ular case of Sagittarius (e.g. Ibata, Gilmore, & Irwin 1994), 
several “ultra-faint” dwarf galaxies exhibit distorted or very 
elongated morphologies (e.g. the Hercules dSph; Roderick et 
al. 2015) however, among the “classical” dSphs, only Carina 
shows clear signs of tidal disturbance in its structural prop¬ 
erties (Munoz et al. 2006; Battaglia et al. 2012; McMonigal 
et al. 2014). Given the difficulties in detecting low-surface 
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brightness tidal features from stellar counts in objects as 
faint and extended on the sky as the MW dSphs, it is still 
not excluded that such non-detections may be due to the 
lack of appropriately deep and spatially extended photomet¬ 
ric data-sets. 

A tumultuous past, marked by strong mass loss due to 
interactions with the host, is expected also from a theoretical 
point of view, in order to explain the present-day properties 
of the observed population of MW satellite galaxies in a 
well-motivated cosmological framework (e.g. Sawala et al. 
2014; Barber et al. 2015), and even more so in those models 
that explain dSphs as the result of “tidally stirred” dwarf 
disc galaxies (e.g. Mayer et al. 2001, Kazantzidis et al. 2011, 
Tomozeiu, Mayer, & Quinn 2015). 

In general, as long as the satellite galaxy retains a bound 
core, the LOS velocity dispersion of the stars in the central 
regions remains a good indicator of the maximum circular 
velocity and bound mass (e.g. Munoz, Majewski, & Johnston 
2008; Penarrubia, Navarro, & McConnachie 2008; Klimen- 
towski et al. 2009; Kazantzidis et al. 2011). At larger pro¬ 
jected radii, however, stars originating in the tidal tails may 
contaminate kinematic samples of individual stars in MW 
dSphs, inflating their measured velocity dispersions and sig- 
nihcantly altering the DM content and distribution inferred 
from an equilibrium dynamical analysis (e.g. Klimentowski 
et al. 2007, Read et al. 2006; Munoz, Majewski, & Johnston 
2008). 

The amount of tidal disturbance experienced is obvi¬ 
ously dependent on the orbital history of the satellite. How¬ 
ever, even if an object has suffered a large degree of mass 
loss, this does not necessarily imply lack of current dynam¬ 
ical equilibrium. Tidally perturbed stars progressively be¬ 
come unbound and are eventually dispersed, with the ob¬ 
ject eventually settling in a new equilibrium configuration 
(Penarrubia et al. 2009). Also, the degree of contamination 
of kinematic data-sets from stars originating in tidal tails 
can vary along the orbit, as a consequence of the varying 
orientation of the inner regions of the tidal tails with re¬ 
spect to the LOS (e.g. Klimentowski et al. 2009). Given the 
difficulty in unambiguously identifying unbound stars lin¬ 
gering close to the main body of the MW dSphs, it is quite 
possible that, depending on the specific orbital history and 
current location along the orbit, the observed kinematics of 
stars is a faithful tracer of their underlying potential in some 
of the MW dSphs and not in others. 

N-body simulations tailored to reproduce the present- 
day position and systemic velocity, as well as the observed 
structural and internal kinematic properties of specific MW 
dSphs appear particularly well suited to explore this issue. 
Such an approach has been for example applied to Carina 
(Munoz et al. 2006) and Leo I (Sohn et al. 2007, Mateo, 
Olszewski, & Walker 2008), for which it was shown that 
modeling these galaxies as tidally disturbed systems initially 
embedded in “mass-follows-light” DM halos, rather than in 
extended halos, yields a good match to the data. The in¬ 
ferred DM content is much lower than what inferred from 
a Jeans analysis in the hypothesis of dynamical equilibrium 
(from 24(M/L)0,\/ in Koch et al. 2007 to ~5 (M/L)©,\/ in 
Sohn et al. 2007, for Leo I; from 80 (M/L) 0 ,v in Walker 
et al. 2007 down to ~40 (M/L)©,^ in Munoz, Majewski, & 
Johnston 2008 for Carina). The orbits that provide the best 
match to the data have small perigalacticons and high eccen¬ 


tricity. While for Carina this is in good agreement with the 
inferences from the measured proper motion (Piatek et al. 
2003), it is likely that Leo I orbit has a larger perigalacticon 
distance, 91±36 kpc (Sohn et al. 2013). Hence the impact of 
tidal interaction with the MW might have been signihcantly 
over-estimated for this galaxy. 

For investigations of the dynamical status of the MW 
dSphs, it is therefore advisable to focus on modeling objects 
placed on observationally motivated orbits. In this work, we 
take under consideration the Fornax dSph. There is consen¬ 
sus that this object is moving on a rather external orbit 
around the MW (Dinescu et al. 2004; Piatek et al. 2007, 
hereafter P07; Mendez et al. 2011), so Fornax is a good can¬ 
didate of dSph that might be close to equilibrium. Our aim 
is to quantify the effects of tides on Fornax’s structure and 
kinematics, by means of A^-body simulations following the 
evolution in the MW of a satellite dwarf galaxy, similar in 
structure to Fornax, on orbits that are consistent with the 
current observational constraints, including Fornax’s proper 
motions. The Fornax dSph is also one of the best studied 
MW dSphs, hence we can then rely on a comprehensive set 
of observational properties against which to compare the 
models. 

We review the relevant set of observables in Section 2. In 
Section 3 we describe the set-up of the Wbody simulations, 
whose results are shown in Section 4. We conclude with a 
summary in Section 5. 


2 THE FORNAX DWARF SPHEROIDAL 

GALAXY: OBSERVED PROPERTIES 

The A-body simulations used in this work have been set 
up in order to reproduce the main observables of the For¬ 
nax dSph. Concerning the structural properties of the stel¬ 
lar component, we focus on reproducing the measured stel¬ 
lar mass and observed surface number density profile, while 
for the internal kinematics we focus on the observed LOS 
velocity-dispersion profile. The simulated object is also re¬ 
quired to have present-day projected position, heliocentric 
distance, systemic radial velocity and systemic proper mo¬ 
tion similar to the observed ones. 

The Fornax dSph appears to have a regular, well- 
behaved 2D stellar distribution, with constant ellipticity 
e ~ 0.3, as revealed by spatially extended photometric data 
reaching down to the base of the red giant branch (RGB; 
e.g. Battaglia et al. 2006; see Bate et al. 2015 for a recent, 
panoramic view over a 25 deg^ region around Fornax based 
on VST/ATLAS data). The only sign of irregularity at large 
spatial scales was claimed by Coleman et al. (2005), who 
detected a shell-like feature outside Fornax’s nominal King 
tidal radius (~ 3.25 kpc) in the distribution of RGB stars 
along the minor axis. The deeper imaging survey by Bate et 
al. (2015) showed that this feature is a mis-identified over¬ 
density of background galaxies. In the inner regions the stel¬ 
lar distribution shows departures from elliptical symmetry 
(Irwin & Hatzidimitriou 1995), which are very likely due 
to the asymmetric spatial distribution of the young main 
sequence stars (see Stetson, Hesser, & Smecker-Hane 1998; 
Battaglia et al. 2006; de Boer et al. 2013). There are also a 
number of small stellar substructures in the spatial distribu¬ 
tion and kinematics of the stars within 0.3 deg (Coleman et 
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al. 2004; Battaglia et al. 2006; Olszewski et al. 2006; de Boer 
et al. 2013), some of which may be related to the accretion 
of a smaller stellar system and others simply to patchy re¬ 
cent star formation. All of the above signs of irregularity do 
not require the occurrence of significant tidal effects. 

Fornax’s distance modulus has been derived from var¬ 
ious techniques - optical photometry of the horizontal 
branch, the tip of the RGB, RR Lyrae stars, red clump (RC) 
stars, as well as from near-infrared photometry of the RC 
and RGB tip (Pietrzyhski et al. 2009 and references therein). 
These are all in good agreement within the typical standard 
and systematic errors, yielding an average value of ~140 
kpc, very close to the value here adopted, i.e. d— 138±8 kpc 
(Saviane, Held, & Bertelli 2000). 

We take as reference the ellipticity e, position angle 
(P.A.) and observed surface number-density profile from 
Battaglia et al. (2006), who used relatively deep ESO/WFI 
V & I photometry extending out to Fornax’s nominal King 
tidal radius to study the structure of this galaxy. The ob¬ 
served surface number density as a function of circularized^ 
radius is best fit by a Sersic (1968) profile 


I(R) = 7o exp 


/_R \ 1/”"' 


( 1 ) 


with Rs = 14.5^ (~ 0.58 kpc, at a distance of 138 kpc) 
and m = 0.71. For these values of Rs and m the projected 
circularized radius containing half of the total number of 
stars is ~ 0.62 kpc. 

The stellar mass of Fornax was estimated by de Boer 
et al. (2012) who performed a star formation history deter¬ 
mination over a wide area using what are so far the deepest 
available color-magnitude diagrams over a large spatial re¬ 
gion for this galaxy. The resulting stellar mass is 4.3x 10^Mq 
within a projected circularized radius 0.66 deg (1.6 kpc). The 
observed LOS velocity dispersion profile was derived by M. 
Breddels from ~2900 individual probable Fornax members 
following the approach outlined in BH13, but binning stars 
in elliptical annuli®. 

A central information to our analysis is the one provided 
by the systemic heliocentric velocity and proper motion of 
Fornax, which determines the orbit of the object around 
the MW and allows us to understand how important tidal 
effects may have been in the evolution of this galaxy. In the 
literature there is good consistency among the heliocentric 
LOS systemic velocity (wsys.h) values derived from different 
spectroscopic data-sets of individual stars, even though the 
samples have widely different sizes and were gathered at 
various observing facilities (e.g. Usys,h = 53± 1.8, 54.1 ±0.5, 
53.3 ± 0.8kms“® from Mateo et al. 1991, Battaglia et al. 
2006, Walker et al. 2006, respectively). 

There are several direct estimates of the systemic proper 


^ In the literature the observed properties of Fornax are often 
reported as functions of the semi-major axis a (also called “ellip¬ 
tical radius”). In this work, in both observations and simulations, 
we always refer to the circularized radius R = \/aE = ay'l — e, 
where b is the semi-minor axis and e is the ellipticity. In particu¬ 
lar, given the observed ellipticity e = 0.3, the circularized Sersic 
radius Rs = 14.5' corresponds to an elliptical radius 17.3'. 

® The LOS velocity dispersion profile published in BH13 differs 
from the one used here, because BH13 binned the data in circular 
annuli. 


Table 1. Values of Fornax’s observational parameters adopted 
in this work. The position angle, P.A., is defined as the angle 
between the North and the projected major axis of the galaxy 
measured anti-clockwise. The systemic velocity Usys^h refers to the 
heliocentric reference frame, while the radial (ur) and tangential 
(nt) velocity components are in the Galactocentric system. The 
quoted radii are circularized projected radii. References: 1) Mateo 
(1998); 2) Battaglia et al. (2006); 3) de Boer et al. (2012); 4) 
Saviane, Held, & Bertelli (2000); 5) Piatek et al. (2007). 
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3 
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4 
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motion of Fornax from astrometric measurements (Dinescu 
et al. 2004; Piatek et al. 2007, hereafter P07; Mendez et 
al. 2011), which are somewhat discrepant (see table 6 of 
Mendez et al. 2011). For example, the current range of most 
likely apogalacticon estimates is very wide, from 152 kpc 
(P07) reaching up to the 965±356 kpc as from the work of 
Mendez et ah, in which case Fornax would be at its first 
pericentric passage. However, all measurements consistently 
suggests that Fornax is unlikely to have passed close to the 
MW centre, as the most likely perigalacticons range from 
118 kpc (P07) to 148 kpc (Mendez et al. 2011). As we wish to 
assess the impact of tidal disturbances, we adopt the orbital 
parameters from P07, which yield the tightest orbit among 
the direct proper motion measurements. This determination 
is also in very good agreement with indirect proper motion 
determinations from stellar redshifts (Walker et al. 2008). 

The values of Fornax’s observational parameters 
adopted in this paper are summarized in Table 1. 


3 SET-UP OF THE Af-BODY SIMULATIONS 


In order to study the effect of tides on the internal struc¬ 
ture and kinematics of Fornax we ran a set Y-body simula¬ 
tions of a Fornax-like object in orbit around the MW. In the 
simulations the MW is represented as a fixed gravitational 
potential: in particular we adopt the Galactic potential of 
Johnston, Spergel, & Hernquist (1995), which consists of a 
Miyamoto & Nagai (1975) disc 


, _ G iVfdisc 


( 2 ) 


a spherical Hernquist (1990) bulge 


0bulge 


G Mbulge 

r + c 


and a logarithmic DM halo 

0halo Uhalo ln(r ± ddalo): 


( 3 ) 

( 4 ) 


where r is the Galactocentric radial spherical coordinate, 
and R and 2 : are the Galactocentric radial and vertical cylin¬ 
drical coordinates. We adopt Mdisc = 10^^ Al©, Msuige ~ 
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Figure 1. Paths of the POTbest (left-hand panels) and P07ecc 
(right-hand panels) orbits integrated in the MW potential for 
12 Gyr up to the present-day position of Fornax. Here x, y and 
2 are Cartesian Galactocentric coordinates in units of kpc (2 is 
orthogonal to the Galactic plane). 


3.4x 10^° Mq, Uhaio = 128kms“^, a = 6.5kpc, h = 0.26kpc, 
c = 0.7 kpc and dhaio = 12 kpc. This Galactic model has 
been previously used in many other studies focused on the 
determination of orbits of globular clusters (Dinescu et al. 
2000; Allen, Moreno, & Pichardo 2008) and by P07 to cal¬ 
culate the orbit of Fornax. 

The phase-space coordinates of the centre of mass of the 
simulated satellite at the beginning of the simulations have 
been derived by taking a test particle with the present-day 
phase-space coordinates of Fornax’s centre and evolving it 
backward in time for 12 Gyr in the adopted Galactic poten¬ 
tial. In order to take into account the uncertainties in the 
estimated proper motions, we considered two different or¬ 
bits: i) the best-fitting orbit inferred by P07 (hereafter, “or¬ 
bit P07best”); ii) the most eccentric (i.e. with the smallest 
perigalacticon) orbit compatible within 3 ct with the observed 
radial velocity and proper motions by P07 (hereafter, “orbit 
P07ecc”). The paths of the two considered orbits are shown 
in Fig. 1: the orbit P07best is an almost polar orbit, with 
perigalacticon rp = 118 kpc and apogalacticon ra = 152 kpc, 
yielding an orbital eccentricity e = (ra —rp)/(ra-|-rp) ~ 0.13; 
the P07ecc orbit has eccentricity e ~ 0.4 and perigalac- 
tic radius rp = 61 kpc. Over the considered 12 Gyr the 
dwarf galaxy goes through three pericentric passages in or¬ 
bit P07best and five pericentric passages in orbit P07ecc. 

In the initial conditions of the simulations Fornax is 
represented with particles as a two-component, spherically 
symmetric system: a stellar component embedded in a DM 
halo. This is an approximation, since Fornax’s stellar com¬ 
ponent is observed to be mildly flattened, with projected 
ellipticity ~ 0.3 (in Section 4.4 we will present a simula¬ 
tion in which we attempt to reproduce also the observed 


ellipticity). The initial structure of the satellite has been 
chosen so that at the end of the simulation (after 12 Gyr) 
the satellite has projected stellar-density and LOS velocity- 
dispersion profiles similar to those observed for Fornax (see 
Section 2). Of course, the problem of matching the final 
snapshot of the simulations with the observations is degen¬ 
erate and requires some assumptions to limit the number of 
free parameters. In particular we assumed that the initial 
velocity distributions of stars and DM in the satellite are 
isotropic, and we also fixed the functional form of the initial 
stellar and DM density profiles. 

The initial stellar component of the dwarf is a spheri¬ 
cal distribution with density profile (Lima Neto, Gerbal, & 
Marquez 1999) 

( ^ 

p* ~ po,* - 

which for u = 1/m and p = 1 — 0.6097i^ -I- 0.05463i^^ is a 
good approximation of the deprojection of the Sersic (1968) 
profile (equation 1) with index m and characteristic radius 
Tc,* = Rs- The density profile of the DM halo is not strongly 
constrained by the available observations: in particular it is 
uncertain whether the central DM density profile is cuspy or 
cored. Though cuspy halos are found in cosmological DM- 
only simulations, cored halos are generally expected in the¬ 
oretical models of dSph formation that account for the ef¬ 
fect of baryons (e.g. Navarro et al. 1996; Mo & Mao 2004; 
Mashchenko et al. 2006; Del Popolo 2009; Di Cintio et al. 
2014; Pontzen & Governato 2014; Nipoti & Binney 2015). 
Moreover, cored halos tend to reproduce the observations 
better than cuspy halos (e.g. Walker & Pefiarrubia 2011; 
Jardel & Gebhardt 2012; Amorisco, Agnello, & Evans 2013), 
though this finding is still debated (see e.g. BH13). Since 
previous works found that the stellar component of simu¬ 
lated galaxies within cuspy DM halos is more resilient to 
tides than in cored DM halos (e.g. Pefiarrubia et al. 2010), 
the assumption of a cored DM profile should enhance the 
pace at which our satellite is affected by tidal disturbances. 
Here we represent the initial DM halo of the dwarf with a 
density profile 

pDM(r) = po.DM (l-\ - ^ e ^ t.Dvi) ^ 

y rc,DM J 

which can be thought of as a NFW-like profile (Navarro, 
Frenk, & White 1995) with scale radius Tc.dm, but with no 
central cusp (din p/d In r = 0 at r = 0) and truncated ex¬ 
ponentially at radii larger than rt,DM. 

In all the simulations we fixed m = 0.71 (the observed 
best-fitting Sersic index; see Section 2), Tc.dm = 3.5 kpc"' 
and rt,DM = 6 kpc, while the initial stellar scale radius rc,*, 



^ The choice of a relatively large value of r„ dm is needed to 
reproduce the flat shape of the LOS velocity-dispersion profile 
observed in Fornax, under our assumption of isotropic stellar ve¬ 
locity distribution. A similar LOS velocity-dispersion profile can 
be obtained with larger values of r,. dmi combined with higher 
values of DM mass, but in practice such haloes would extend far 
beyond the tidal limit and lose a large fraction of their particles 
at the beginning of the simulation. For this reason, and to limit 
the number of free parameters, we decided to keep the value fo 
^■c.DM fixed. 
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total stellar mass = 47r r^p*(r)dr and total DM mass 
Mdm = dTT ff° r^pDM(r)dr are left as free parameters. 

The values of Tc,*, M* and Mom are tuned in each 
simulation to reproduce in the last snapshot the projected 
surface-density and LOS velocity-dispersion profiles of For¬ 
nax, and its observationally estimated stellar mass. For this 
purpose, we used the following iterative procedure. 

(i) We first adopt guess values of rc,*, and Mom that 
give good fits to the present-day projected surface-density 
and LOS velocity-dispersion profiles of Fornax. 

(ii) When a simulation with these values of the param¬ 

eters is run, its final snapshot is projected on the plane of 
the sky and used to measure the stellar mass within a pro¬ 
jected radius of 1.6 kpc the half-light radius (f?™) 

and the LOS velocity dispersion within the half-light radius 
('^LOS.hl)- 

(iii) New values of rc,*, M« and Mom (here indicated with 
primes) are then calculated as follows: 

r-;* = rc,* f 

Ml = M, {M:%/Ml%ff 

Mom = Mom (o'LOS,hi/o'LOS,hi)i 

where = 4.3 x = 0.6 kpc and cr^^^s.hi = 

12kms“^ are the observational estimates for Fornax (see 
Section 2). 

Steps (ii) and (iii) are repeated until convergence (de¬ 
fined by the condition that the three free parameters change 
by less than 10%). For the simulations presented here, the 
above scheme converged after 2-3 iterations. It must be 
noted that our iterative procedure is not as robust as other 
more complex algorithms developed for this purpose (see 
e.g. Ural et al. 2015) and could potentially fail to converge if 
significant changes in the structure of the satellite occurred. 
However, given the relatively wide orbit of Fornax, in all our 
simulations the three parameters change by small amounts 
during the entire satellite evolution (see Section 4), so we 
believe that the simple procedure adopted here provides re¬ 
liable results in the present context. 

The N-body simulations have been performed with 
the collisionless code fvfps (Londrillo et al. 2003; Nipoti 
et al. 2003), which has been specifically modified to in¬ 
clude the fixed external potential of the MW (equations 2- 
4). We adopted the following values of the code parame¬ 
ters: minimum value of the opening parameter Omin = 0.5, 
softening parameter e = 10 pc, and (constant) time-step 
At = 3.2 X 10® yr. The initial equilibrium two-component 
satellite system is generated assuming isotropic velocity dis¬ 
tribution: given the density distribution of each component 
and the corresponding total gravitational potential, the par¬ 
ticles are assigned phase-space coordinates x, y, z, Vx, Vy, 
Vz (centred in the dwarf galaxy’s centre of mass) using the 
ergodic distribution function obtained with Eddington’s in¬ 
version (Binney & Tremaine 2008). The parameters of the 
simulations, including the number of stellar (A«) and DM 
(Ndm) particles, are given in Table 2. 


4 RESULTS 

4.1 Evolution of the satellite’s stellar and 
dark-matter components 

We present here the results of three simulations (P07best, 
P07ecc and P07eccsh) differing in either the satellite orbit 
or the initial stellar and DM density distributions (see Ta¬ 
ble 2). Simulations P07best and P07ecc have the same ini¬ 
tial DM and stellar distribution and differ only in the orbital 
parameters (e = 0.13 and e = 0.4, respectively). In simula¬ 
tion P07eccsh the satellite is placed on the more eccentric 
orbit (as in simulation P07ecc), but the stellar component 
starts with a more centrally concentrated spatial distribu¬ 
tion than in simulation P07ecc, and a DM mass higher by 
a factor 1.4. In all these simulations the total stellar mass 
is M* = 5 X IO^Mq. The simulations are run for 12 Gyr: 
the final snapshot corresponds to the present day and can 
be compared with the observational data of Fornax. 

Figures 2-4 show the spatial distribution of DM and 
stellar particles after 2, 7 and 12 Gyr for the simulations 
P07best, P07ecc and P07eccsh- In all cases the DM compo¬ 
nent develops tidal tails: the fraction of particles contained 
in the tidal tails is clearly higher in the simulations on the 
more eccentric orbit. Tidal tails are instead absent in the 
stellar component in all these simulations. This can be quali¬ 
tatively understood considering that along the P07best orbit 
the minimum (pericentric) tidal radius, calculated according 
to the prescriptions of Allen, Moreno, & Pichardo (2006, 
see their appendix A) and assuming a mass for Fornax of 
3 X 10®Mq (see Table 2), is rt,min = 13.16 kpc, which is 
large compared to the characteristic size of the stellar com¬ 
ponent of Fornax (in the initial conditions of the simulations 
about 99% of the stellar mass is contained within 2kpc). 
The more eccentric orbit P07ecc has a minimum tidal radius 
Ft,min = 7.36 kpc, which is significantly smaller than that of 
orbit P07best, but still larger than the radius containing 
most of the stellar mass. So, during the satellite orbit, stars 
are always contained within the instantaneous tidal radius 
and are not expected to be significantly affected by tides. 
Taken at face value, the above estimates of the tidal radii 
are relative large also with respect to the characteristic size 
of the DM distribution (in the initial conditions about 90% 
of the DM mass is contained within r = 7.5kpc), so the 
above simple analysis would predict moderate tidal effects 
also for the halo. In fact, the A-body simulations show that 
the DM halo develops evident tidal tails and lose a signifi¬ 
cant fraction of its mass. The stellar component is not tidally 
stripped, but, at least on the more eccentric orbit, its inter¬ 
nal structure and kinematics evolve due to the variation of 
the gravitational potential induced by the tidal stripping of 
the halo. 

Figures 5-7 show the intrinsic angle-averaged stellar and 
DM density profiles of the satellite in simulations P07best, 
p07ecc and p07ecCah at t = 2, 7 and 12 Gyr, compared to 
the initial ones. Figures 8-10 show, for the same simulations, 
the initial {t — 0) and final (t = 12 Gyr) projected number- 
density and LOS velocity-dispersion profiles compared to 
the observations. The final projected properties of the satel¬ 
lite are computed by assuming that the system is observed 
from the Sun. Given that we consider an axisymmetric MW 
potential, we have the freedom to fix the position of the Sun 
a posteriori in a circle of radius 8.5 kpc lying in the Galactic 
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Table 2. Initial parameters of the simulations, rc,*: scale radius of the stellar distribution; M^: stellar mass; Mdm- DM mass; Mtot- 
total mass; N^: number of stellar particles; number of DM particles; Ntot'- total number of particles; e: orbit eccentricity. Notes. 

For simulation POTbestgat, in which the stellar component is initially flattened, rc,* refers to the spherical distribution used in the set-up 
of the initial conditions (see Section 4.4). In simulations POTbestjviFL a-nd POTeccjviFL there is no distinction between stars and DM, so 
rc,* is the scale radius of the total mass distribution. 


Name 

rc,* 

kpc 

M* 

Mq 

Mdm 

Mq 

^tot 

Me 

Af* 

Ndm 

Wot 

e 

P07best 

0.581 

5 X 10'^ 

3 X 10® 

3.05 X 10® 

51200 

1536000 

1587200 

0.13 

P07ecc 

0.581 

5 X 10'^ 

3 X 10® 

3.05 X 10® 

51200 

1536000 

1587200 

0.40 

P07ecCsii 

0.465 

5 X 10’^ 

4.2 X 10® 

4.25 X 10® 

51200 

2150400 

2201600 

0.40 

P07best]viFL 

0.581 

- 

- 

1.5 X 10® 

- 

- 

51200 

0.13 

P07eccMFL 

0.581 

- 

- 

1.5 X 10® 

- 

- 

51200 

0.40 

P07bestflat 

0.581 

5 X 10'^ 

2.6 X 10® 

2.65 X 10® 

51200 

1331200 

1382400 

0.13 


plane and centred in the Galactic centre. The Sun position is 
chosen in each simulation to obtain the best possible match 
with the heliocentric centre-of-mass position and systemic 
velocity of Fornax which are in all the simulations presented 
here reproduced to within 5%. 

The final properties of the dwarf galaxies in the sim¬ 
ulations are summarized in Table 3. In simulation P07best 
neither the DM nor the stellar distributions evolve substan¬ 
tially. At the end of the simulation 97.5% of the DM parti¬ 
cles and virtually all the stellar particles are bound to the 
satellite. At the end of the simulation both the projected 
number-density and LOS velocity-dispersion profiles of the 
stars are very similar to the initial ones and compare well 
to the observations. A convenient reference radius to ana¬ 
lyze the properties of the satellite is 1.6 kpc, which contains 
about 90% of the stellar mass in the initial conditions, be¬ 
sides corresponding approximately to the outermost point 
of the observed LOS velocity dispersion profile of Fornax, 
within which de Boer et al. (2012) estimated Fornax’s stel¬ 
lar mass. From the start to the end of the simulation, the 
DM mass within 1.6 kpc changes little, from 2.68 x 10® 
to 2.64 X 10® Mq. The final stellar mass within 1.6 kpc is 
M*,i.6 — 4.6 X lO^M©, consistent within 7% with the esti¬ 
mate of de Boer et al. (2012). 

A satellite with the same initial conditions, placed onto 
the more eccentric orbit (simulation P07ecc), experiences a 
significant evolution due to tidal forces: the surface-density 
profile of the stars becomes more extended and the average 
LOS velocity dispersion becomes lower (see Fig. 9). These 
changes are mainly due to the loss of mass from the DM 
halo, whose density profile is substantially altered over 12 
Gyr of evolution (see Fig. 6). At the end of the simulation 
the satellite retains only 50% of the DM particles (the DM 
mass within 1.6 kpc is 59% of the initial value). Because of 
the consequent variation of the gravitational potential, the 
satellite reaches an equilibrium configuration characterized 
by a larger scale length and lower velocity dispersion (the 
stellar component retains 92% of its mass within 1.6 kpc). 
A similar behaviour was noticed in previous works study¬ 
ing the effects of tides on dwarf galaxies (e.g. Penarrubia, 
Navarro, & McConnachie 2008). Overall, the final snapshot 
of simulation P07best is not a good representation of the ob¬ 
servations: the stellar particle surface-density profile of the 
simulated galaxy is too extended with respect to the ob¬ 
served one, and the LOS velocity-dispersion profile is too 
low (see Fig. 9). 


In order to reconcile with the observations the final 
properties of our simulated galaxy when placed on the more 
eccentric orbit, the initial conditions had to be modified, 
starting with a more concentrated stellar density profile and 
increasing the initial DM mass (model P07ecCsh). As shown 
in Fig. 10, the final snapshot of this simulation is in good 
agreement with the observations. In this case, the fraction of 
bound DM particles at the end of the simulation is 71% (the 
final DM mass within 1.6 kpc is 79% of the initial value). As 
in simulations P07best and P07ecc, virtually all the stellar 
particles are bound to the object after 12 Gyr of evolution. 
The hnal stellar mass within 1.6 kpc is more than 98% of the 
initial value; the object expands slightly, with an associated 
small decrease in the central LOS velocity dispersion. 

The results of simulations P07best and P07ecCsh show 
that the present-day observed structure and kinematics of 
Fornax can be reproduced by a system whose stellar com¬ 
ponent varies only little during its 12 Gyr journey in the 
MW on observationally-motivated orbits. In particular, the 
simulated systems are found to be quasi-stationary in the 
last Gyr of evolution. This suggests that Fornax can be rea¬ 
sonably well modeled as an isolated system in equilibrium. 
The hnal fraction of unbound stars in the simulations is neg¬ 
ligible, so the spectroscopic samples of Fornax stars should 
not be contaminated signihcantly by tidally stripped stars. 
Interestingly, baryons and DM appear to contribute simi¬ 
larly to the mass budget in the central regions of Fornax, 
suggesting that this galaxy may not be “DM dominated at 
all radii”, as believed to be the case for most MW dSphs. 

4.2 Jeans analysis of the final snapshots 

As a further test of the dynamical status of Fornax, we car¬ 
ried out an equilibrium mass modeling analysis of the h- 
nal snapshots of simulations P07best and P07ecCsh, which 
were shown to match well the observed properties of Fornax. 
Should the stellar component of the simulated systems be 
out of dynamical equilibrium, then we would expect to hnd 
a discrepancy between the DM mass determined from the 
equilibrium analysis and the “true” DM mass as measured 
directly from the snapshot. 

The simulated stellar component is our tracer of the 
gravitational potential: given the LOS velocity-dispersion 
and surface number-density profiles measured in the final 
snapshots of the simulations (Figs. 8 and 10), we wish to 
hnd the total mass distribution using the spherical Jeans 
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equations (Binney & Mamon 1982) under the assumption 
that the velocity distribution of the tracers is isotropic. We 
model the total mass density profile as the sum of the stellar 
and DM density profiles. The tracers’ density and mass pro¬ 
files are obtained by deprojecting (assuming spherical sym¬ 
metry) the Sersic profile (equation 1) with parameters as in 
Table 1 and stellar mass 5 x 10^ M©. For the DM density 
profile we assume a simple model of the form 


pDM{r) = pO,DM 



(7) 


(Dehnen 1993), which yields the analytic mass profile 


MuM{r) = Mdm,( 


rc,DM + r 


where Mdm.oo = drrpo.DMr® dm/3 is the total DM mass. 

We determine the best-fitting Mdm.oo and Tc.dm via a 
minimization. Though these two parameters are quite 
degenerate, we find that the DM mass enclosed within 1.6 
kpc (Mdm.i.s) is not significantly affected by this degen¬ 
eracy (Mdm.oo is anyway a ill-constrained quantity due to 
the limited spatial extent of the tracer kinematic data with 
respect to the extent of the DM halo). We obtain a remark¬ 
ably good agreement between Mdm,i.6 yielded by the equi¬ 
librium Jeans analysis and the values measured in the final 
snapshots of simulations P07best and P07ecCsh, with a dif¬ 
ference of at most 7%. This lends support to the fact that, 
at the present-day position of Fornax and along orbits com¬ 
patible with literature proper motion measurements, it is 
reasonable to assume that tidal effects are not significant 
for DM determinations of this satellite galaxy within the 
region probed so far by the kinematic data. 


4.3 Mass-follows-light models 

So far we have assumed that the initial satellite’s DM halo is 
extended. In principle, an alternative possibility is that, in 
the absence of an extended DM halo, tidal forces are effective 
in producing a fiat and high LOS velocity-dispersion profile 
as observed (Aaronson 1983). For instance, such a model 
appears to provide a good match to the observed structure 
and kinematics of the Carina dSph (Munoz, Majewski, & 
Johnston 2008). To explore this idea for the case of Fornax, 
we ran two additional simulations in which we make the ex¬ 
treme assumption that “mass follows light”. In practice, in 
these simulations we do not distinguish between stars and 
DM: the satellite is represented initially as a one-component 
isotropic system of total mass Mtot = 1.5 x 10®M© with the 
same initial density profile as the stellar component of sim¬ 
ulations P07best and P07ecc. With this choice the central 
LOS velocity dispersion is similar to the value measured in 
Fornax. However, as expected, the initial velocity-dispersion 
profile is much steeper than observed. We considered both 
the P07best orbit (simulation POJbestMFu) and the more 
eccentric orbit P07ecc (simulation P07eccMFL). The result¬ 
ing spatial distribution of particles after 2, 7 and 12 Gyr are 
shown in Figs. 11 and 12; the corresponding intrinsic den¬ 
sity profiles are shown in Fig. 13, and the surface-density 
and LOS velocity-dispersion profiles in Fig. 14. 

In simulation P07bestMFL a small fraction (~2%) of the 
initial mass appears to be decoupled from the main body of 


the galaxy, thus forming a hint of tidal tails. At the end of 
the evolution both the projected number-density and LOS 
velocity-dispersion profiles of the stars are very similar to 
the initial ones. Therefore, the LOS velocity dispersion pro¬ 
file is still significantly decreasing with radius, being just 
somewhat higher than the initial profile at i? > 0.6 deg, 
with velocity dispersion ~ 6kms“^. This is at odds with 
observations, in which the velocity dispersion profile is flat 
with values « 12 km/s up to J? ~ 0.7 deg (see Fig. 14, left- 
hand panels). In simulation P07eccMFL the system develops 
a clearer pair of tidal tails. However, also in this case the 
satellite loses a small fraction of the initial satellite mass 
budget (~3%), so the final density and velocity dispersion 
profiles are very similar to those of simulation P07bestMFL. 
High values of the final LOS velocity dispersion are found 
in simulation POLoccmfl only at 77 > 1 deg, far beyond the 
region covered by kinematic observations (see Fig. 14, right- 
hand panels). We note that in the orbits considered here 
the tidal tails developed by the satellite are never aligned 
along the LOS (a configuration that could enhance the LOS 
velocity dispersion at smaller projected radii; Klimentowski 
et al. 2007). 

These results suggest that, even in the absence of an 
extended DM halo, neither the large amplitude of Fornax’s 
velocity-dispersion profile, nor its approximately flat shape, 
are caused by lack of dynamical equilibrium due to tidal dis¬ 
turbances from the MW. In other words, Fornax’s kinematic 
samples are not expected to be signihcantly contaminated 
by tidally stripped stars. 

4.4 A flattened model 

In the simulations described in the previous sections we have 
assumed that the satellite is initially spherical. In fact, we 
mentioned in Section 2 that the observed stellar component 
of Fornax is flattened, with observed ellipticity e — 0.3. The 
simplest interpretation of this deviation from circular sym¬ 
metry is that Fornax’s stellar component is intrinsically ax- 
isymmetric. Given the available observational information, 
the inclination of the system with respect to the LOS is 
unknown and it is unclear how much rotation or velocity 
anisotropy contribute to the observed flattening. A more 
realistic simulation of the evolution of Fornax in the MW 
should start from a non-spherically symmetric satellite and 
match, in the final snapshot, the observed ellipticity, to¬ 
gether with the other observables. Clearly this is not an 
easy task, because intrinsic shape, inclination, rotation and 
anisotropy add several degrees of freedom. In the spirit of 
showing that the observed ellipticity can be reproduced suf- 
hciently well with a simple model and that the results de¬ 
scribed above are not specihc to spherically symmetric mod¬ 
els, here we present, as an example, the result of an addi¬ 
tional simulation named P07bestfiat, which follows the orbit 
P07best (see Table 2). 

In simulation P07bestflat the satellite is initially repre¬ 
sented by an oblate (but non-rotating) stellar component 
and by an almost-spherical DM halo. In practice, the initial 
conditions are constructed as follows. We start by generating 
a two-component spherical, isotropic system with stellar and 
DM density prohles given by equations (5) and (6), respec¬ 
tively, with the same parameters as in simulation P07best, 
but Mdm = 2.6 X 10® Mq. Each particle is characterized by 
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coordinates x, y, z, Vx, Vy, Vz centred in the dwarf galaxy’s 
centre of mass. The phase-space coordinates of the stellar® 
particles are modified as follows: x, y, Vx and Vy are multi¬ 
plied by a coefficient a, and « and Vz by a coefficient b < a. 
In practice we use a = and b = so the final 

minor-to-major intrinsic axis ratio is q = b/a. Specifically, 
we adopt q = 0.58, so the (constant) intrinsic ellipticity is 
Ci = 1 — g = 0.42. Such a system is not in equilibrium, but 
we verified empirically that it is actually a quasi-equilibrium 
configuration. We followed the evolution of such a system in 
isolation with an Wbody simulation, finding that an equilib¬ 
rium configuration, still with ellipticity ei ~ 0.41, is reached 
after 1 Gyr of evolution. The DM distribution slightly ad¬ 
justs itself because of the presence of the flattened stellar 
component, but we find that after 1 Gyr the halo is still 
very close to spherical symmetry. We take as initial {t = 0) 
satellite in simulation P07bestflat this A^-body system as it is 
after 1 Gyr of evolution in isolation. The initial phase-space 
coordinates of the centre of mass of the dwarf are the same 
as in simulation P07best and the orientation of the system’s 
axes is empirically chosen so that at the end of the simula¬ 
tion the system is seen almost edge-on when observed from 
the Sun. The iterative algorithm described in Section 3 has 
been used to find the optimal parameters at the beginning 
of the simulation, including the initial intrinsic ellipticity Ci, 
which has been updated at each iteration according to the 
formula 

f / obs / sim\0.5 

/« ) , 

where = 0.3 and are the projected ellipticities cal¬ 
culated on observed data and in the last snapshot of the 
simulation (we took as initial guess ei = 0.3). The initial 
orientation of the system is such that the equatorial plane is 
parallel to the LOS from the Galactic centre to the observed 
position of Fornax. During the simulation the vector normal 
to the satellite’s equatorial plane changes its direction, os¬ 
cillating by several degrees, but in the last snapshot the 
satellite is still seen from the Sun position as approximately 
edge-on. The final properties of simulation P07bestflat are 
reported in Table 3. 

Figure 15 shows the distribution of stellar particles in 
the tangent plane, centred at the coordinates of the object 
centre of mass, as seen from the Sun. The projected elliptic¬ 
ity of this simulated object is e ~ 0.3, in perfect agreement 
with the value measured for Fornax. As it can be noted 
from Fig. 16 both the resulting surface number density and 
LOS velocity profile are also in this case in relatively good 
agreement with Fornax’s observed properties. The DM and 
stellar mass loss is similar to the P07best spherical case: less 
than 3% of stars and ~8% of DM are lost within a radius 
of 1.6 kpc. In Fig. 16 we also overlay the profiles of the sim¬ 
ulated stellar component at the end of the evolution, as it 
would appear if seen exactly edge-on: in this case the ob¬ 
ject’s ellipticity would be e ~ 0.42, but the surface-density 
and LOS velocity dispersion profiles would not be much dif¬ 
ferent from those found projecting along the actual LOS in 
the simulation. 


® The phase-space coordinates of the DM particles are left un¬ 
changed. 


Table 3. Final properties of the dwarf galaxy in the two- 
component simulations. M* i.g: stellar mass within an angle- 
averaged radius r = 1.6kpc; A7 dm,i. 6: DM mass within r = 
1.6 kpc; final-to-initial stellar mass ratio within r = 1.6 kpc; 

/oM.l.e: final-to-initial DM mass ratio within r = 1.6kpc. 


Name 

A7*,i,6 
10'^ Mq 

A7dm,1.6 
10® Mq 

/*,1.6 

/dm,1.6 

P07best 

4.62 

2.64 

0.991 

0.985 

P07ecc 

4.30 

1.58 

0.923 

0.590 

P07ecCsii 

4.70 

2.95 

0.975 

0.787 

P07bestflat 

4.53 

2.40 

0.971 

0.917 


5 DISCUSSION AND CONCLUSIONS 

In this work we simulated the dynamical evolution over a 12 
Gyr orbit around the MW of a satellite, whose final struc¬ 
ture and kinematics resemble those of the Fornax dSph. The 
results shown in Section 4 indicate that even assuming a 
relatively eccentric orbit (still within 3a the best-fit sys¬ 
temic velocities estimated by P07) the effects produced by 
the Galactic tidal field are relatively small, and the final sys¬ 
tem is close to equilibrium in its own gravitational potential. 
This suggests that the assumption of dynamical equilibrium 
can be safely used to estimate the mass (and DM content) 
of Fornax and that the existing kinematic samples are not 
expected to be contaminated by tidally stripped stars. 

Our result is robust in the sense that our models were 
built so to maximize the possible tidal effects, because we 
followed Fornax orbit for as long as 12 Gyr in a “fully grown” 
MW, and we considered the orbits with the smallest peri- 
centre still consistent with direct measurements of Fornax 
proper motion. According to hierarchical models of galaxy 
formation, a system like the MW is expected to have sub¬ 
stantially increased its DM mass over 12 Gyr and, of course, 
also the MW stellar components have not been in place as 
they are now for such a long time span. Furthermore, For¬ 
nax may have become a MW satellite more recently than 
12 Gyr ago. For example, Rocha, Peter, & Bullock (2012), 
based on the results of the Via Lactea II simulation, propose 
an infall time for Fornax of 5-9 Gyr ago. Another indication 
that the orbital evolution of Fornax in the MW was likely 
shorter than 12 Gyr can be found in its very extended star 
formation history, with star formation having occurred un¬ 
til ~50-100 Myr ago. It is still debated on what timescales 
ram-pressure stripping can remove the gaseous component of 
the MW dSphs progenitors; nonetheless, the observational 
finding that Local Group gas-rich star-forming dwarfs are 
mostly located beyond the estimated MW or M31 virial ra¬ 
dius (e.g. van den Bergh 2000, Irwin et al. 2007) suggests 
that objects with recent star formation activity must have 
spent a considerable fraction of their lifetime far from the 
large Local Group spirals. Overall, we expect that in a more 
realistic set-up the effect of tides would be milder than what 
examined here. 

As we describe the MW as a fixed gravitational po¬ 
tential, we neglect the effect of dynamical friction on the 
satellite’s orbit. As Fornax orbits in relatively low-density 
regions of the MW, dynamical friction is not expected to be 
very important. In any case, the effect of dynamical friction 
is to bring the satellite closer to the host’s centre: if dynam- 
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ical friction were accounted for, the backward in time recon¬ 
structed orbit of Fornax would be more extended than the 
one considered here. Also in this sense, our simulations with 
fixed MW potential tend to slightly overestimate the tidal 
effects. Moreover, a more external orbit and consequently 
weaker tidal effects would be obtained by assuming a less 
massive model for the MW, as recently suggested by Gib¬ 
bons, Belokurov & Evans (2014) on the basis of the shape 
and extent of the Sagittarius stream. 

Based on the results of the present work, we are re¬ 
assured on the fact that dynamical mass estimates of For¬ 
nax derived from equilibrium dynamical models are reliable. 
Still, as well known, dynamical models suffer from the mass- 
anisotropy degeneracy, so that different combinations of DM 
density and stellar-anisotropy profiles can reproduce equally 
well the observations. Moreover, given Fornax’s observed el- 
lipticity, additional degrees of freedom derive from the un¬ 
certainties on the intrinsic shape and rotational support of 
this dSph (see Hayashi & Chiba 2012 and Jardel & Geb- 
hardt 2012, who studied axisymmetric equilibrium models 
of Fornax). However, while the inner DM distribution is hard 
to determine, the mass within 1.6 kpc (approximately the 
outermost point of the observed LOS velocity dispersion pro¬ 
file of Fornax) is believed to be more robustly constrained 
(BH13). Therefore, in the limits of the considered approx¬ 
imations (initial spherical symmetry and isotropic velocity 
distribution, with the exception of simulation P07bestflat), 
our simulations can be also used to estimate Fornax’s DM 
mass over the scales probed by the existing kinematic data¬ 
sets. The simulations that successfully reproduce Fornax’s 
observables have a dark-to-luminous mass ratio within 1.6 
kpc in the range 5 — 6 (see Table 3). Measuring the same 
ratio within 3 kpc yields values in the range 16—18. We con¬ 
firm that Fornax is amongst the least DM dominated dSphs, 
and we suggest that baryons and DM may contribute simi¬ 
larly to the mass budget in Fornax’s central regions. In the 
near future we plan to extend the present analysis to other 
dSphs for which the orbital parameters are sufficiently well 
constrained. 
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Figure 2. Distribution of DM (black) and stellar (red) particles of the satellite galaxy as seen in the x-y, y-z, and x-z planes (from left 
to right) for simulation POTbest, after 2, 7 and 12 Gyr of evolution (from top to bottom). Here x, y and 2 are Cartesian Galactocentric 
coordinates in units of kpc (z is orthogonal to the Galactic plane). For convenience we plot one every thirty DM particles. 
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Figure 3. Same as Fig. 2, but for simulation POTecc. 
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Figure 4. Same as Fig. 2, but for simulation POTecCgh 
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Figure 5. Angle-averaged density profiles of DM (black) and stellar (red) particles for three snapshots, at 2, 7 and 12 Gyr, for the 
simulation POTbest (solid lines) as compared to the initial conditions (dashed lines). 


DM Initiol 
DM t= N Gyr 
Stors Initiol 
Stors t= N Gyr 



Figure 6. Same as Fig. 5, but for simulation P07ecc. 
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Figure 7. Same as Fig. 5, but for simulation P07eccsh- 
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Figure 8. Surface number-density (normalized to the central 
value; top panel) and LOS velocity-dispersion (bottom panel) 
profiles of the stellar component of the POTbest model. The or¬ 
ange and black lines indicate the initial conditions (t=0) and the 
final properties (t = 12 Gyr), respectively. The observed sur¬ 
face number-density profile from Battaglia et al. (2006) and LOS 
velocity-dispersion profile from BH13 are shown as black asterisks 
with error bars. The grey band in the bottom panel refers to the 
LOS velocity-dispersion profile from the 12 Gyr snapshot when 
using a similar number of stellar particles and binning as done by 
BH13 for the observations (3000 particles, 250 per bin). 




Figure 9. Same as Fig. 8, but for the P07ecc model. 






Figure 10. Same as Fig. 8, but for the POTecCgh model. 
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Figure 11. Same as Fig. 2, but for simulation POTbestMFL- Here mass follows light, so we do not distinguish between stellar and DM 
particles. 
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POybestypL 



POyecc^FL 



Figure 13. Angle-averaged density profiles for four snapshots, at 0, 2, 7 and 12 Gyr, for the POTbestMFL (left) and P07eccMFL (right) 
models. Here mass follows light, so we do not distinguish between stars and DM. 






Figure 14. As in Fig. 8, but for the P07bestMFL (left) and P07eccMFL (right) model. 
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Figure 15. Spatial distribution of stellar particles in simulation 
P07bestflat projected in the tangent plane as seen from the Sun. 
The ellipse has the ellipticity, position angle and nominal King 
tidal radius as in measured by B06 (see Table 1). 




Figure 16. As in Fig. 8, but for the P07bestfiat model. The black 
curves show the profiles as seen from the Sun. The orange and 
blue lines show the profiles of the initial and final snapshots, as 
seen edge-on. 
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